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L.  summary  and  Introduction 

y  A  theory  aad  computer  program  are  presented  for  computing  the  transonic 
fan  duct  and  free  jet  flow  exhausting  over  the  centerbody  of  a  turbofan 
engine  and  merging  with  an  exterior  stream  flowing  over  a  nacelle .^The" 
physical  problem  under  consideration  Is  sketched  in  Figure  1.  For 
simplicity,  thejprimary  jet  is  treated  as  a  solid  circular  body.  The 
exterior  stream  and  fan  duct  flows  are  both  inviscid  and  are  represented 
by  appropriate  differential  equations  based  on  nonlinear  transonic  small 
disturbance  theory.  The  position  of  the  free  jet  boundary  is  found  as 
part  of  the  solution  by  imposing  the  free  slip  boundary  condition  of  con¬ 
tinuity  of  pressure  and  slope  at  the  linearized  boundary  interface.  The 
Mach  number  of  the  expanded  jet  at  downstream  infinity  where  the  static 
pressure  is  ambient  is  specified,  and  this,  in  turn,  establishes  the 
total  pressure  of  the  jet  flow  by  means  of  the  isentropic  relation.  The 
total  jet  mass  flux  is  not  a  free  parameter  but  is  determined  as  a  part 
of  the  solution  by  the  choice  of  the  downstream  jet  Mach  number  and 
application  of  the  Kutta  condition  at  the  trailing  edge  of  the  fan  cowl. 
The  mass  flux  at  the  fan  plane  (see  Fig.  1)  is  assumed  to  be  uniformily 
distributed.  Also  required  is  the  distribution  of  axial  velocity  along 
a  transverse  plane  exterior  to  the  nacelle  and  upstream  of  the  fan  exit 
(see  Fig.  1),  and  the  Mach  number  of  the  undisturbed  external  flow  at 
infinity. 

Numerical  solutions  are  obtained  by  finite  difference  algorithms  based  on 
the  original  developments  of  Murman,  Cole,  and  Krupp  (Refs.  1-9). 

Figure  2b  shows  a  typical  example  of  computed  streamline  patterns,  Jet 
boundary  slope  and  sonic  lines.  For  this  particular  example,  there 
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FIGURE  1  ILLUSTRATION  OF  THE  MODEL  PROBLEM 


Ot  4100  7740  OHIO.*/ 71  JIB-047 


FIGURE  2  PRESSURE  DISTRIBUTION,  SONIC  LINES,  AND  STREAMLINE 
PATTERN  FOR  EXAMPUZ  FLOW 
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occurred  an  imbedded  region  of  subsonic  flow  adjoining  the  centerbody 
downstream  of  the  fan  cowl  exit  plane.  Figure  2a  shows  the  corresponding 
pressure  distributions  on  the  centerbody,  on  the  interior  wall  of  the 
fan  duct  and  along  the  free  streamline  dividing  the  jet  and  the  exterior 
flow.  No  comparisons  have  yet  been  made  with  experimental  data. 

This  report  presents  a  complete  derivation  of  the  theory,  a  description 
of  the  finite  difference  procedures  employed,  and  several  examples  of 
calculated  flows.  A  description  of  the  computer  program  is  given  in  the 
Appendices  together  with  instructions  for  its  use. 


2.  Derivation  of  the  Transonic  Small  Perturbation  Differential  Equations  for 


the  Flow  in  the  Duct  and  in  the  External  Stream 


Since  the  entropy  change  across  a  weak  shock  is  third  order  in  the 
pressure  jump,  we  shall  assume  irrotational  flow  and  introduce  a  velocity 
potential  which  satisfies  the  following  differential  equation: 

k  (°}-$*)  *fr,  r  Ca-'~&  ’  -  2  £  +<£$n  /r-0!1) 

where  x  and  are  the  dimensionless  axial  and  radial  coordinates  based  on 
duct  length,  and  a  is  the  local  velocity  of  sound  given  by 


with  denoting  the  velocity  at  infinity  downstream.  We  now  introduce 

a  perturbation  velocity  potential  if  with  a  small  parameter  £  in  the 


£  =  Li,(x  +e<f) 


Since  small  disturbances  are  propagated  normal  to  the  principal  flow 
direction  in  a  nearly  sonic  free  stream,  we  introduce  a  scaling  factor  in 
the  r^  direction.  This  is  expressed  by 

r  *  tr, 

where  X  snd  fr  will  be  determined  appropriately  in  the  following  develop¬ 
ment. 
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Boundary  conditions  on  the  body  are  linearized.  Let  r^(x)  =  d/^+SR^{x) 
define  the  boundary  of  the  nacelle  adjoining  the  free  stream.  Here  is 
a  thickness  ratio  associated  with  either  the  nacelle  or  the  centerbody.  At 
r^  ■  d/2,  the  boundary  condition  that  the  flow  be  tangential  to  the  nacelle 
wall  becomes 

h  //*  * 

Substituting  the  perturbation  velocity  potential  yields 

n<fr/v  !3) 

It  is  convenient  to  let 

x  -  <r/t  oo 

and  the  boundary  condition  takes  the  familiar  linearized  form 

-  7?JoO  r=  r rt  s'tJ./z  (5) 

The  velocity  of  sound  to  the  first  order  in  6  becomes 

>/n%  ■+ 0(t‘)  (6) 

2  2 

Substituting  the  preceding  relation  for  a  into  the  differential 

equation  (1)  and  neglecting  terms  of  the  order  of  S  yield 

■»  £%  £  Vrr  +  VrJ=  a4J<7) 

We  consider  the  limit  as  approaches  unity  while  the  perturbation 

parameter  6  approaches  zero.  Ihus,  we  define  the  limit 

)/ert£  K 

Letting 
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■*  =  cm! 


(8) 


renders  the  j*  derivatives  of  the  same  order  as  the  x  derivative  and  yields 
the  well  known  transonic  small  perturbation  differential  equation  in  the 
form 


O  -cr+i)<fl  yxx  •+  (r<{r)rlr  =  o 


(9) 


Equations  (8)  and  (k)  enable  us  to  determine  f  and  “C  in  terms  of 
the  free  stream  Mach  number  and  thickness  ratio  £  ;  i.e., 

.*/»  r//3 


and 


(1C) 


*/» 


(11) 


Also  for  X,  we  have 


m  /  r  v/s 

KzO-M*)/*  (~\m 

Note  that  in  the  limiting  process  we  retained  some  non-linearity  in 
the  differential  equation.  Equation  (9)  is  elliptic  when 

<fx  <  /r/rr+O 

and  hyperbolic  when 

Yk  >  */cr+') 

The  quantityjK/(y+  t )  is"  the  sonic  velocity. 

We  now  apply  a  similar  approach  to  the  duct  flow.  Let  M.  be  the  Mach 

J 

number  at  infinity  downstream  where  the  duct  flow  and  outside  stream  have 
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the  same  static  pressure.  Then  for  the  perturbation  velocity  potential  we 
define 


j*  =  Uj  (*■  ■*  «,•  y  ) 


The  velocity  of  sound  in  the  duct  is  given  by  a  relation  similar  to  Eq.  (6) 


0?!  Uj  *  *  /"* -Cr-t)6j*fx  +0(6/  J  (13) 

Substituting  Eqs.  (12)  and  (13)  into  the  differential  equation  (l)  leads  to 
the  equation 

”  o(*J 

i  i  J  J  (1U) 


We  shall  retain  the  same  scaling  factor  T  as  in  the  outside  stream. 


=■  /m*  (16^ 

The  transonic  parameter  for  the  duct  flow  is 

nt  «  \*  -  O-rtj)  )  U7) 

and  the  differential  equation  for  flow  inside  the  duct  becomes 

-«■  -  o  <l8> 

For  the  duct  walls  defined  by  r^  ■  d/2  +  J" R^(x) ,  the  linearized  boundary 
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conditions  take  the  form 


or 

Yr  ~  <?0< 

The  transonic  small  perturbation  differential  Eqs.  (9)  and  (18)  are 
non-linear  and  practical  analytical  solutions  are  difficult  to  obtain.  In 
the  sequel,  the  equations  are  expressed  in  difference  fora  and  numerical  re¬ 
laxation  techniques  are  applied  to  find  the  solutions  for  specific  flow  con¬ 
figurations  . 


3.  The  Pressure  Coefficient  for  the  Transonic  Small  Perturbation  Theo: 


Pt  4100  7740  OWIft.l/TI _  JIB -04  7 


Following  a  similar  procedure  for  the  duct  flow  we  obtain  for  the  pressure 
coefficient, 

’  (p  -  -Z6j  <fx  - 

Since  the  pressures  of  the  two  streams  balance  at  infinity,  ?.  =  P 

W 

and  the  condition  that  the  pressure  balance  on  the  boundary  streamline 

between  the  two  flows  is  easily  seen  to  be 

_  Af>0) 

x  x 

where  the  superscripts  ce>  and  j  denote  the  limits  of  freestream  and  duct 
velocity  perturbation  potentials,  respectively,  as  one  approaches  the 
boundary  streamline. 


4.  Pate  of  Mass  Flow  in  the  Stream  Direction 

The  magnitudes  of  at  the  fan  plane  evolves  as  a  part  of  the  solution 
and  is  subject  to  periodic  updating  during  the  course  of  the  relaxa¬ 
tion  procedure.  For  this  purpose,  we  derive  a  formula  for  the  mass 
flow  rate  which  is  required  for  updating  ^  as  the  computation  pro¬ 
gresses.  Neglecting  the  cross  component  of  velocity,  we  obtain  the 
following  relation  for  the  local  density  in  the  jet 


* 


Expanding  and  retaining  terms  to  the  order  of  £•  yields 


« 


+ 


of*/) 

(21) 


This  may  be  written  in  the  following  form  by  regrouping  the  second  order 


terms: 


REV  SYM 


MO,  D6  -41078 

l6 

pin  -- '  -  m!*j  i*  -(r  ->j  r> 

\ ,  .  .  .  +*;*>;-o&/a 

Since  My6yf/^-  “V  -  —  Cy  /(  - Ofcjthe  third  term  may  be  neglected. 


The  rate  of  mass  flow  per  unit  area  then  is  given  by 

/>/*  A.;  4*  /  + 

/'  ■+  •*  °(*f) 


We  note  that  this  may  be  written 


ff* Ifrt  '•  '  <2u> 

"/  + 

The  last  term  is  seen  to  be  of  higher  order  than  the  preceding  terms.  The 

reason  for  retaining  second  order  terms  in  becomes  apparent  when  we 


introduce 


ejO-Mf)  -  6jMyKl  *•  «  */r. 


/*  *  *  /  f  <>(*;  ) 


into  Eq.  (24)  and  obtain 


p£  /fi°j  =  -(r+dfjtl 


Note  that  the  mass  flow  is  a  maximum  when  the  flow  is  sonic,  i.e.,  at 

<ft  *  K,  /er+< ) 

Solving  Equation  ( 25)  for  yields 


ft  -  f  ^  l/fr-n) 

(fix.  -fius)//>Jus- 


where  a 
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With  the  mass  flow  rate  a  computed  from  an  intermediate  solution,  the 
value  for  ip  at  the  fan  plane  is  updated  in  terms  of  the  values  of  <p* 
from  Equation  (27),  keeping  the  total  mass  flux  in  the  jet  fixed.  .A  rela¬ 
tion  similar  to  3q.  (27)  for  the  exterior  stream  is  easily  written  down 

by  replacing  K.  and  e  by  X  and  e ,  respectively. 

J 

5.  The  Physical  Problem  to  be  Solved  and  the  Development  of  the  Bounds 


Value  Problem 

We  are  interested  in  finding  the  pressure  distribution  over  the  back  part 
of  a  nacelle  and  on  the  duct  centerbody  over  which  the  air  from  the  fan 
flows.  The  problem  is  idealized  to  the  flow  field  illustrated  in  Figure  1 
with  a  solid  body  replacing  the  central  jet  flow.  Since  we  linearize  the 
boundary  conditions  on  the  solid  boundaries  and  on  the  free  streamline 
boundary  between  the  duct  flow  and  the  outside  stream,  the  idealized  flow 
is  expressed  as  a  boundary  value  problem  for  the  perturbation  velocity 
potential.  This  is  presented  in  Figure  3» 


The  boundary  value  problem  appears  to  be  well  posed,  since  the  value  of 
the  potential  or  its  normal  derivative  is  prescribed  on  the  boundary.  How¬ 
ever,  the  application  of  a  Kutta  condition  at  the  fan  cowl  trailing' edge  is 
an  additional  restriction  which  requires  that  the  magnitude  of  at  the 
fan  plane  not  be  arbitrary,  but  must  be  determined  from  the  solution  of 
the  problem.  This  value  is  found  from  computing  the  mass  flow  in  the  free 
jet  far  downstream  and  determining  the  value  of  <P _  at  the  fan  plane  which 
satisfies  conservation  of  mass  in  the  duct  and  jet.  The  mass  flow  from  the 
duct  calculated  far  downstream  is  not  a  sensitive  function  of  the  initial 
guess  for  <p*  at  the  fan  plane  (which  must  be  input)  but  is  dependent  ulti¬ 
mately  upon  the  choice  of  undisturbed  free  jet  Mach  number  because 
of  the  Kutta  condition.  The  input  fan 
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plane  boundary  conditions  are  corrected  from  calculations  of  the  con¬ 
servation  of  mass  as  the  solution  progresses.  The  frees tream  Mach  number 

and  the  free  jet  Mach  number.  My  are  two  parameters  which  must  be  fixed 
in  the  calculations. 


The  differential  equations  in  the  two  regions  are  expressed  in  dif¬ 
ference  form  and  the  flow  field  computed  by  means  of  relaxation  techniques. 


6.  The  Difference  Equation  for  the  Interior  Points  of  the  Flow 


To  apply  relaxation  techniques  to  the  boundary  value  problem  in  Fig.  3, 
we  divide  the  region  into  a  rectangular  mesh  by  lines  parallel  to  the  x  and 
r  axes.  The  spacing  is  generally  not  uniform  as  greater  accuracy  is  achieved 
in  minimum  computing  time  by  using  a  finer  spacing  near  boundaries  for  the  r 
variable  and  finer  spacing  near  rapidly  accelerating  or  decelerating  regions 
of  the  flow  in  the  x  variable.  Let  the  subscripts  i,  j  denote  the  values  of 
the  variables  at  typical  interior  points  given  by  the  coordinates  x^,ry 
The  second  derivative  in  the  x  direction  at  x  s  x^  and  r  •  is  expressed 
in  difference  form  for  subsonic  flow  by  using  the  values  at  points  on  the 
both  sides  of  the  point  i,  J  and  is  easily  seen  to  be 

=  t  -in.  -  ~ ^ x 

*,,,  -x,  ,  <  V*.  *«->) 

,  (28) 

r  -  te,  Y,j 


c-i  -  //<**, 

~  J /  C  ~  X i )  i  -  Xj-i  ) 


v  ' 1 J  (■  \  '  Vx.-  -  *•> } 
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The  notation  as  used  here  is  suggestive  of  the  notation  used  in  the  Fortran 
statements  of  the  computer  programs  in  order  to  make  it  easy  for  the  user 
to  understand  the  computer  programs  from  this  discussion. 


For  the  derivative  with  respect  to  r,  we  use  a  similar  central  difference 
formula,  for  both  subsonic  and  supersonic  flow. 


(r<f,v/r=  { q-xaiy.-j a  ~V  -  n-ti.cf.-i 

rs*.  -1  -• c r,  -  \  1  ri(rj+rrs-. 

~  2<tj  ^<J -/  -Zfo-j  +  6j)¥lS  *  Zbi  fij-'  (3C 


where 


*j  ’C-'A  A;-  oj+(- '/..Vo  ~C--. ) 


(3D 


and  we  introduce  the  convenient  notation 


(32) 


To  complete  the  difference  fora  of  the  differential  equation,  we 
require  a  representation  of  the  first  derivative  with  respect  to  x.  For 
points  equally  spaced  a  xlX  apart,  the  difference  form 

<fu,j  -  Avj  )  /*** 

approximates  ^  to  the  second  order  in  <XX  at  the  central  point  i,j. 
To  obtain  a  second  order  fora  for  ^  when  the  points  are  not  equally 
spaced,  we  write 


"O.D6-41D78 

e*at  21 
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?X  -  C/ 2  (%Hj  ~%'i*  * *  j  j  )  (33) 

and  determine  C/ and  d.  f  {  appropriately  by  Taylors  itnci  expansions 
about  ij.  Let  h^  *  xi+1  -x^  and  h2  =  x.^  15:1611  expanding  about  the 


point  i , J  yields 


■+  ~  ^Vxx/i  *~)  ~  ft. 

Setting  the  coefficient  of  ^  equal  to  1  and  the  coefficient  of  the 
term  equal  to  zero  yields  two  linear  equations  to  solve  for  and  V 


The  result  is 


C,  i  *  A  <W  'f rrirlu  -x v  \  ' c'' 

fV,/  r>-)(jri><  'V'->  (35) 

4,- 

For  inside  the  duct  we  define  the  difference  form 

l/-.  r  O-/*lj)/r/0-A*Z)  -(*  +  •)?*  (36) 


and  for  the  outside  stream,  we  use 


2it.j  c  /T  -  (37)  — 
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By  applying  Sqa.  (28),  (30),  and  (36)  to  the  differential  equation  (18) 
when  the  flow  in  the  duct  is  subsonic,  we  obtain  the  following  difference 


equation 


*  vrj  ( V,+,S  -e>  *>)  * 

-  a  t »/ Y.y.,  ~‘*i  +-*j)  %■  ■*  ® 


In  the  relaxation  process,  we  solve  for  the  values  of  along  a  column 
x  *  in  the  flow.  It  is  therefore  acre  convenient  to  write  Eq.  (38)  as 

-  "  * rj  (ci  ^r*tj  +  ^ 

In  this  form,  with  v^  regarded  as  fixed,  a  set  of  linear  equations  for  the 
values  of  along  x  »  xt  results.  The  matrix  of  coefficients  is  tri¬ 

diagonal  and  is  easy  to  solve.  The  non-linear  equation  is  solved  by  itera¬ 
tion.  The  quantity  v^  is  defined  by  the  values  of  <f  from  the  previously 
iteration  and  then  is  found  from  solving  Eq.  (39).  The  calculations 
are  repeated  with  v^  determined  by  the  recently  obtained  values  of  .  The 
iteration  is  continued  until  the  required  accuracy  is  achieved. 

When  the  flow  is  supersonic,  i.e.,  when  v^.  <  0,  we  use  a  backward 

difference  formula  in  the  axial  variable  x  instead  of  the  central  difference 
formula  described  above. 


Hence, 

[*,  -fr+i)  £  J  =  V,„j ((,-1  f!-j  " 
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With  Eq.  (^0)  substituted  for  the  first  term  in  Eq.  (38),  we  obtain 
the  following  equation  equivalent  to  Eq.  (39)  for  supersonic  flow; 

“  +  (e>-,  j  ?-/ ¥*-*y  ) 


The  points  at  which  v  <,  0  but  v^  ^  >  0  are  parabolic .  When  this  occurs , 

Eq.  (4l)  is  modified  by  setting  =0. 


For  the  interior  points  of  the  outside  stream,  similar  equations  result 
but  with  u.  and  u.  replacing  v  and  v  ,  respectively. 


7.  The  Application  of  Wall  Boundary  Conditions  to  the  Difference  Equation  at 


Mesh  Points  Adjacent  to  Boundaries 


The  application  of  the  boundary  conditions  near  a  solid  wall  is  more 
easily  made  when  the  mesh  points  are  not  included  on  the  boundary  but  a  half 
mesh  point  width  away.  Consider  first  the  lower  boundary.  Figure  4  shows 
the  points  i,  j  *  1  and  2  and  the  boundary  point  where  the  boundary  con¬ 
dition  is  applied.  The  r  derivative  becomes 

1  fr  (P  }  -  -  (r  V  -  tAWflZlL*  ^  f 

r(rTr>f  “  rt(r0-r%)iri>‘r  r,  -  r%  J  (**) 


where  for  convenience  we  choose  rQ  ■  r^  -  (r^  -r^)  *  2rfe  -r^.  From  Fig.  3, 

the  introduction  of  the  boundary  condition  yields 


REVSYM 


MO.D6-41078 
•>*oi  24 


AO  IM«D 


REV  SYM 


FIGURE  4  COLUMN  MESH  POINTS  ILLUSTRATING  THE  APPLICATION  OF 
THE  LOWER  BOUNDARY  CONDITION 


z 


FIGURE  5  COLUM!  MESH  POINTS  ILLUSTRATING  THE  APPLICATION 

OF  THE  UPPER  BOUNDARY  CONDITIONS  AND  THE  CONTINUITY 
ACROSS  THE  FREE  STREAMLINE 
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where 


■‘rirf^r  • 


_  2 r^M 

lb,  Iftx  -  Y,V  ) 


1.  -  <rbMj  />-l”2trx-r.) 

The  difference  form  of  the  differential  equation  then  becomes 

vij I*. , w 

Rewriting  in  tridiagonal  form  for  the  ith  column  relaxation  yields 

"A  +  V,*  -  (u5) 

-*ii  (CfYf+n 

When  the  point  is  supersonic,  we  have  instead, 

c  ^V_,  /  £•*-•  ~  i ?r3i  -  (46) 

K-»  /  ^V- /  ^ t-i  r-z  &{ 

For  the  upper  boundary,  situated  at  r  *  r „,  halfway  between  the  r. 

and  r.  ,  we  obtain  a  similar  relation  of  the  r  derivative  (see  Fig.  j). 
Jm  +  1 

A.-'w  I  I 

-  v  ri+rKr*>  (u7) 

’'m  '** 
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where  «x  -  “  rf<t )  ^  . 

The  equation  equivalent  to  Eq.  (*+5)  then  becomes 

G-jJfij-i  ci  %'l  r 

When  v.  .  ^  0,  then  Eq.  (49)  becomes 
m 

+  ( -*-J  * V 


(48) 


(49) 


(50) 


For  the  point  r  =  r .  the  application  of  the  boundary  condition 

^m 

in  Fig.  2  by  comparison  with  Eq.  (45)  is  seen  to  yield 


+  bu<  ‘'•it i 

-  ~  ui  J*i  Ce>  i  + 

to 


(c-i 


where  **  *  6^*  -  rf  -J)  .  (52) 

A  formula  for  the  hyperbolic  case  analogous  to  Eq.  (46)  for  the  outside 


stream  may  easily  be  written  down. 


8.  Application  of  the  Prescribed  Inlet  Flow  Condition  to  the  Difference  Equation 

Along  the  line  x  *  0  in  the  duct  as  well  as  in  the  outside  stream,  we 
prescribe  the  values  of  ^  at  the  points  of  the  mesh  r  -  r^,  j  *  1  to 
jMr.  Let  s  f(r)  at  x  *  0  and  define  f^  *  f(rj).  For  ft  at  x  *  x^, 
we  follow  a  procedure  similar  to  that  used  for  the  r  boundary  conditions. 
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For  convenience  we  introduce  the  potential 


Y-  c  ftx  -rr+>)Y 

Then  the  x  derivative  term  of  the  differential  equation  becomes 

LK  V'XJf  /(*+() 

-  -fafx  !*■(** o 


Expressing  this  quantity  in  difference  form  leads  to 


where  the  superscripts  3/2  and  0  denote  the  value  of  V*  at  x  =  (x^  +  x^J/2 
and  x  =  0,  respectively.  We  define  xQ  and  we  obtain,  after  factoring 

Eq.  <5*0, 

-y*  i/c*t+x.) 


Since  <fx  =  f  , 


we  have 


X  "  'x 


ft -(*+•)  £• 


Replacing  the  variable  ^  yields 


define  fy]  - 


To  make  the  notation  consistent  with  the  nomenclature  for  general  values  of 
x,  we  write 

u.j  fe. /?*;-'/’./)  -*4,4/1 

(58) 

where  cu  =  1/2  (Xg  -x^  and  d^  =  l/2.  The  differential  equation  at  the 
point  x  *  and  r  =  r^  takes  the  following  difference  form 


JL  x?-x* 


-ill 


-  D 


Expressed  in  tridiagonal  form,  this  becomes 

t)  tC<  “/;)  *«;  (60) 

-  **< /  (4  tj  -  c,  'ft.i ) 

2  2 

where  c^  =  l/(x  ^  -  x  ^)  and  *  l/(Xg  +  x^).  Similar  equations  may 

readily  be  written  for  the  duct  equation  with  v^  defined  by  Eq.  (56)  with 

k  replaced  by  (1  -  M^)k/(1  -  M^).  The  relation  for  j  =  1  and  j  =  Jm  +  1 

is  easily  found  by  setting  *  0  and  adding  the  appropriate  boundary  term 

to  the  right  hand  side.  Similarly  for  j  *  J  ,  b.  is  set  equal  to  zero  and 
t  m  jm 

the  ldded  to  the  right  hand  side. 


For  subsonic  points  in  the  flow  field,  over -relaxation  is  employed  by 
the  program  in  the  manner  of  Murman  (8).  However,  it  was  found  that  over¬ 
relaxation  at  the  inlet  interfered  with  the  satisfaction  of  the  ^  boundary 
condition  at  x  *  0.  This  was  remedied  by  varying  the  relaxation  parameter 
w  continuously  from  unity  at  x  -  0  to  about  1.3  to  1.7  at  x  ■  X(10). 


-O.D6-41Q71 
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Conditions  for  the  Interface  Between  the  Free 


Stream  and  the  Duct  Flow 


The  boundary  conditions  of  continuity  of  slope  and  pressure  across  the 
free  streamline  boundary  will  be  applied  at  the  line 


r  -  1  - 


for  x  x  or  i^iQ.  We  have  shown  that  the  relations  between  the  stream- 
10 

line  slope  and  .the  r  derivative  for  both  the  interior  and  exterior  flows  are 

VT  •  < 

Continuity  of  slope  across  the  streamline  then  yields 


rc*:0-  "jc' 


For  the  continuity  of  pressure  we  have 


«)  . 


Since  this  relation  must  hold  at  all  points  of  the  free  streamline,  we  can 
integrate  from  the  trailing  edge  point  to  the  point  x^ .  We  obtain 

</7x  - )  -  f  (%  )  -  -  tM(x ) 

Letting  r  *  7'  we  restate  our  boundary  condition  as 

'  *•  *0 

rf)~i r, )  -  *  (63) 


we  restate  our  boundary  condition  as 


REV  SYM 


Dl  4100  7740  OHIO. «/T  I 


The  quantity  Wijf  is  a  constant  for  the  entire  free  streamline  and  is  computed 
from  the  values  at  the  trailing  edge  at  each  iteration. 


We  shall  satisfy  the  boundary  conditions  on  the  free  streamline  at 

the  centerline  point  between  r  =  r .  and  r  =  r.  , .  To  simplify  the  re- 

m  in 

lations  we  shall  further  assume  that  the  lines  r=r.  ,  to  r  =  r.  „ 

jm-l  jm+2 

are  equally  spaced  a  distance  h.  Since  is  defined  only  at  the  mesh 
points  we  apply  Taylor's  expansions  to  find  at  the  midpoint  r  =  rf  = 
rjm+l/2‘  Thus >  for  ^  on  both  sides  of  the  free  streamline,  we  have 

Kb. 

(6U) 

-  V,-;  H-lfJ  ♦ 

-  ri'?„  3 

where  we  have  dropped  the  superscripts  on  the  right  hand  side  since  the  sub¬ 
script  adequately  defines  the  values.  In  the  preceding  two  equations 
we  replace  and  derivatives  by  differences.  As  in  Eqs.  (33)  through 

(35)  for  the  x  derivative,  we  develop  a  second  order  difference  for  V  and 


obtain 


—  r,,-  _,')/**' 

r  jr*r*  *• 

Similarly  for  the  second  derivatives,  we  obtain 

„  I  „*>{.-*  — 

Trr,i}*.  L  r  J- 

=  -  V  ^ 

*rr(r_  r 
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Substitution  of  Eqs.  (65)  and  (66)  into  Eqs.  (64)  and  solving  for  <f 
and  y^yields  the  following  two  equations  for  the  values  of  on  the 


two  sides  of  the  streamline 

*•) 


'8  4  9  f,u>  -  **  Y* 

9KJ8  -Y.  jJ*  *  1  •>  Yr°/S 


/» 


(67) 


,  f-gp)  es 

Two  more  equations  for  Jr  and  can  be  obtained  from  the 

difference  forms  of  the  partial  differential  equations  at  the  mesh  points 

i,j  and  i,j  +1.  The  r  derivative  for  r  *  r.  takes  a  form  similar  to  Eqs. 
mm  Jm 

(42)  and  (47);  i.e., 


CJJ 


Jr,  J~ 


(68) 


From  the  definition  of  a  and  b  we  recognize  that  Eq.  (68)  may  be  written 

a) 


\  (rfjr  ]  *  i  u.4jjfr  -  7, •/_-,)]  (69) 

Similarly  for  r  =  r  +, ,  we  find  that 
•'m 

|  Crfr)r  |  *  8  f (70) 

fs  f' 

r  0,4.1 

Substituting  Eqs.  (69)  and  (70)  into  the  difference  form  of  the  differential 
equations  yields,  for  v^  and  u^ 
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(71) 

u"(c<f*  -  -£.?•■ +*( ,7.  .  )->?«* {••  -V,v,,)-° 

Uij+)lW)  +  \  *'*-')*'  Jff  Jm  JJX  J* 

^"0  (fl0‘) 

3r  later  convenience,  we  solve  Eqs.  (71)  for  I  ^  and  7  and  rearrange 


the  terms  to  obtain 


hbyJrJ>T  +A*‘f,>-  *  BL 

+A'*'*‘  *  3' 


where 


A>* -(eivijjr  A 

B,  *  J,*') 

*»*  ''’Le'  +*j~ 

6t  --  -  ,W.+  **>  ■> 

When  the  flow  on  either  side  is  supersonic  C^ij  or*})  )  J 

then  Eq.  (72)  holds  but  with  A^,  B^  and  and  B2  defined  by 

■  ^*1-1  ~ 

0,  z“u/-‘j+ter-t '^l“' ^»“*/>/ 1  (?u , 

5  ^ f  I’w  ^ ^ 

With  the  boundary  conditions  of  Eqs.  (62)  and  (63),  the  quantities 

*fj<**^  *f  ^  nay  te  eliminated  from  Eqs.  (67)  and  (72). 
Using  Eqs.  (62)  in  the  form 
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%  w  •• 


(75) 


yields  from  Eq.  (72) 


*  Y*7J 

O.. 


(76) 


When  we  subtract  Eqs.  (67)  we  obtain  the  single  equation 

*  7,  9  %r%rih(rrJ* 


(77) 


where  Af  is  a  known  constant  which  for  the  final  converged  solution  must 
agree  with  the  jump  at  the  nacelle  trailing  edge  and  at  infinity  downstream. 
Eliminating  the  r  derivatives  in  Eq.  (77)  by  Eqs.  (72)  leads  to 

£l  ■+' 

WMf  ■ 

(78) 


3  7,'jV. ,  -  f  9  ♦  5  Vu.)  7, • +  f 9  -  M,  /<*;_.  , 

-  c<-*3 ■*  6*/b<‘J  -8A<f  =  ° 


Equations  (76)  and  (78)  are  the  equations  for  j  ;  JB  and  j  =  j^+1 
required  to  complete  the  set  of  equations  for  the  column  relaxation  procedure. 
In  this  form,  the  matrix  is  no  longer  tridiagonal  and  more  sophisticated 
methods  for  solution  would  be  required.  However,  the  tridiagonal  matrix  is 
restored  by  first  eliminating  from  Eq.  (76)  and  (78)  and  then 

eliminating  f  .  After  simplification,  there  results 
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(79) 


(80) 


r«;bjp 

+  +  ,,  +  3^,1]$  <fVj 

+  V  9 .My  +  M;  byA  ,1  <f  ,-Jm4  , 

-r -  SK$ b;Jb;£?.o 

+ 1 9  / n  l  a-j^,  -  a  ,i  3  Mi  a.jM’ m  *  r s  ^  b;„  y j}  , 

+^9,,^ p  -  f5/li  V  bjJti'B,  SM&fj&'f*  o 

Equations  (79)  and  (80)  are  based  on  a  second  order  expansion  of  ^ 
about  the  streamline  points.  A  simpler  first  order  theory  can  also  be 
derived.  In  place  of  Eqs.  (6U)  we  write 

V(“=  ~  £  Ly^  0'O 

With  given  by  Eqs.  (65),  we  have 


Similarly, 


V01:  (3?.-;  -7,  -  ,)/« 


'/W1 


REV  SYM 

OOSTAMO 

*w.D6-din7ft 

M«  35 

9- 


Ol  4  100  7740  OHIO. «/7  t 


I 


Substracting  these  last  two  equations  yields 

This  relation  is  simpler  than  Eq.  (67)  since  it  does  not  contain  the  first 
derivative  terms.  The  continuity  of  slope  leads  to  Eq.  (76)  for  the  other 
equations  required  to  complete  the  system  of  column  equations.  The  tri¬ 
diagonal  form  of  the  system  is  preserved  by  first  eliminating  - 

and  then  . *  ,  from  the  two  equations.  This  yields  the  following 

two  equations  in  place  of  Eqs.  (79)  and  (80) : 

f  <  V/  ^  \  bL+>  f>is  • 

—  (Mj  b’  T5,  )  “  ZMj  bj  bj 

yvy  ***  ^  ^ 

-^V'% 

10.  Application  of  the  Method  of  Relaxation  to  the  Solution  of  the  Boundary 
Value  Problem  in  Figure  3. 

In  the  foregoing  analysis,  we  have  formulated  the  difference  equations 
corresponding  to  the  partial  differential  equations  of  Eqs.  (9)  and  (18) 
for  points  in  the  interior  of  duct  and  outside  stream,  for  mesh  points  near 
the  solid  boundary,  on  both  sides  of  the  common  streamline  of  the  two  flows, 
and  near  the  inlet  regions.  We  have  expressed  these  equations  in  a  form  for 
solving  the  values  of  along  a  column  of  fixed  values  of  i.  These  equation! 
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are  non-linear  in  since  the  quantities  u^  and  also  include  the 
term  *ffj  .  When  u^  and  v^  are  assumed  to  be  fixed  and  determined  by 


the  starting  or  previously  iterated  values  of  ,  then  the  set  of 

equations  for  the  1  column  is  linear  in  ■  The  matrix  is  tridiagonal  in 
form,  and  is  especially  easy  to  solve.  Each  column  is  iterated  until  con¬ 
secutive  approximations  of  ‘fjj  are  within  a  specified  accuracy.  For  a 
more  complete  description  of  the  numerical  procedure  the  reader  is  referred 
to  the  papers  by  Murman,  Cole  and  Krupp  in  the  references at  the  end  of 


The  required  relations  for  computing  the  flow  field  have  been  presented 
in  sufficient  detail  in  the  preceding  analysis  to  set  up  the  procedures  for 
the  numerical  evaluation  of  the  flow  field.  In  Appendix  I,  the  results  of 
the  foregoing  analysis  will  be  summarized  in  a  form  used  directly  in  the  re¬ 
laxation  procedure.  In  this  way,  the  reader  will  be  able  to  follow  more 
easily  the  computer  program  coded  to  perform  these  calculations. 


11.  Calculation  of  Streamline  Pattern  and  Duct  Mass  Flow 

The  application  of  the  Kutta  condition  at  the  nacelle  trailing  edge 
requires  that  ^  at  the  fan  plane  not  be  arbitrary  but  be  determined  as 
part  of  the  solution  for  a  specified  value  of  far  downstream.  This  Is 
found  by  requiring  that  the  total  mass  flow  at  the  duct  inlet  be  equal  to 
the  calculated  far  downstream  value.  For  convenience  we  prescribe  a  uniform 
mass  flux  distribution  at  the  fan  plane  and  aasvae  that  the  streasOines 
far  downstream  are  also  parallel.  The  mass  flow  rate  is  then  given  in 
terms  of  ^  by  Eq.  (25 ) . 

To  evaluate  the  mass  flow  rate  in  the  free  jet  far  downstream  we  must 
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where  AQ  is  the  boundary  area  of  the  duet  between  r^  and  rj  , 

2q.  (27)  yields  the  appropriate  value  of  <fK  to  assign  to  the  inlet  for 
the  next  iterations.  The  inlet  value  of  ^  is  updated  from  time  to  time 
until  the  relaxation  converges  to  suitable  accuracy. 

12.  Examples  of  Calculated  Flows 

To  demonstrate  the  operation  of  the  program  TEA-343,  the  configuration 
described  by  the  coordinates  in  Tables  1,  2,  and  3  was  used.  The  maximum 
radius  of  the  centerbody  wall  and  the  end  of  the  fan  cowl  and  outer  duct 
wall  all  occur  at  the  axial  position  x  =  1  in  the  dimensionless  variables. 
The  thickness  ratio  d  used  in  the  parameters  K  and  e  was  chosen  as  the  in¬ 
crease  in  centerbody  radius  from  x  =  0  to  the  maximum  radius  at  x  =  1 
divided  by  the  section  length.  For  the  data  in  Tables  1,  2,  and  3,  this 
value  is  S=  0.15.  For  convenience,  the  boundary  value  at  x  =  0  in  the 
duct  was  chosen  as  uniform  over  the  fan  plane.  The  boundary  values  along 
x  =  0  in  the  outside  stream  were  chosen  to  vary  with  r  like  l^r^+l)^2 
which  is  the  variation  for  the  incompressible  solution  fro  a  source  and 
3ink  at  the  axial  positions  x  =  +  1  in  a  uniform  stream  given  by  Milne- 

Thorason  (10)  on  page  486.  With  on  the  fan  cowl  at  x  =  0  given  by  C  = 

^  P 

-0.06  the  flow  was  computed  for  a  Mach  number  of  the  exterior  stream  of 

M„  =0.9  and  a  duct  far  downstream  Mach  number  M.  =  0.9.  The  duct  Mach 
oo  *  J 

number  Mj  was  increased  from  Mj  =  0.9  to  1.25  by  increments  of  .05.  Each 
step  required  200  to  300  iterations  and  complete  convergence  was  not 
necessarily  achieved  at  each  step.  It  was  found  that  the  iterations  failed 
to  converge  if  the  first  estimate  of  at  the  inlet  of  the  duct  was  too 
far  from  the  correct  value.  A  first  estimate  of Is  found  by  satisfying 
the  conservation  of  mass  flow  through  the  fan  plane  under  the  assumption 
of  a  uniform  flow  at  Mach 
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TABLE  1 


COORDINATES  AND  SLOPE  OF  NACELLE 


I 

X 

R 
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COORDINATES  AND  SLOPE  OF  UPPER  FAN  CONTOUR 
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COORDINATES  AND  SLOPE  OF  LOWER  FAN  CONTOUR' 
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number  Mj  through  the  fan  exit  plane  at  the  trailing  edge  of 
the  fan  cowl  (Fig.l),  and  the  program  corrects  this  value 
during  the  computations. 

The  distributions  of  pressure  coefficient  on  the  centerbody 
and  on  the  duct  flow  outer  boundary  are  presented  in  Fig.  6  for 
Mj  =  1.2.  The  coefficients  are  based  on  the  exterior  free 
stream  dynamic  pressure.  The  velocity  is  seen  to  be  highest 
and  the  pressure  coefficient  lowest  at  the  point  of  maximum 
centerbody  radius  x  *  1.  The  streamline  pattern  of  this  flow 
is  shown  in  Figure  7.  The  duct  flow  is  accelerated  to  supersonic 
velocities  leading  to  the  sharply  curved  sonic  line  terminating 
close  to  the  fan  cowl  trailing  edge.  Note  that  the  flow  near 
the  centerbody  is  rapidly  decelerated  by  a  shock  resulting  in  a 
small  local  region  of  subsonic  flow. 

Figure  8  shows  the  pressure  distribution  on  the  centerbody 
and  duct  flow  free  streamline  for  a  duct  Mach  number  of  Mqo  =  1.2 
For  this  value  of  the  Mach  number,  the  duct  flow  remains  super¬ 
sonic  downstream  of  the  curved  sonic  line  as  seen  in  Figure  9. 

The  flow  decelerates  from  the  maximum  value  at  x  »  1  but  does 
not  become  subsonic.  The  streamline  pattern  is  very  similar  to 
that  for  Mqq  **  1.20. 

The  centerbody  of  Figure  9  was  extended  downstream  as  shown 
in  the  streamline  pattern  of  Figure  10  (See  Table  3) .  This 
corresponds  more  closely  to  the  actual  flow  conditions  and  the 
resulting  pressure  distribution  shown  in  Figure  11  is  similar  to 
that  of  Figure  8  in  the  vicinity  of  the  fan  exit  plane. 
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DISTRIBUTION  OF  PRESSURE  COEFFICIENTS  ON 
THE  CENTERBODY  AND  FREE  STREAM  BOUNDARY. 
PRESSURE  COEFFICIENT  IS  BASED  ON  DYNAMIC 
PRESSURE  OF  FREE  STREAM. 
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FIGURE  9  STREAMLINE  PATTERN  IN  THE  FAN  DUCT  WITH  SONIC  LINES  FOR  FLOW  CONDITIONS 
OF  FIGURE  8 
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FIGURE  11  DISTRIBUTION  OF  PRESSURE  COEFFICIENTS  ON  THE  CENTERBOBY  AND  FREE  STREAM 
BOUNDARY  OF  FIGURE  10 


APPENDIX  I:  SUMMARY  OF  THE  EQUATIONS  USED  IN  THE 
COMPUTER  PROGRAM  TRANSDUCT 

The  equations  to  be  solved  at  each  column  i  have  been  expressed  in  a 
form  which  yields  a  tridiagonal  N  x  N  matrix  for  the  coefficients  of  the 
variables,  where  N  is  the  number  of  values  of  (p  in  a  single  column.  The 
elements  of  the  matrix  are  described  by  three  N  dimensional  vectors.  The 
coefficient  of  the  diagonal  term  ^  in  the  Jtb  equation  is  the  com¬ 
ponent  of  the  vector  denoted  by 

DIAG(J) 

in  the  program.  The  vector  of  coefficients  of  ^  to  the  right  of  the 

diagonal  is  designated  as 

SUPER(J) 

while  the  vector  of  coefficients  to  the  left  of  the  diagonal  Cp  ^  is 
designated  by 

SUB(J) 

1.  Coefficients  Based  on  the  Mesh  Points 

When  the  distributions  of  the  x^  and  r^  are  established,  the  co¬ 
efficients  depending  upon  these  quantities  can  be  determined  once  and  for 
all  at  the  ouset  of  computations.  For  the  x  coordinate,  we  define 

cL  - 1/(4  -  4) 

ci  *  V^iei  ’  xi>  <Xi+l  •  Xi-1> 

*1  -  VU2  -  V 

di  *  V<xi  *  Xi-1>  (xi+l  •  Xi-1} 

cu-  1/2(x2  -  -  (*t  -  ximl) 

*11  1/2,^  -  (xi  +1  -  xt)  dt 


(Si) 


U/l‘»IHO  0»lt  001*  IQ 


For  tbe  r  coordinate,  we  obtained 


*j  '  (ri  *  r3-i)/2  ra  (ra*i  '  rj-i>  (ra  '  rj-i>  3  >  1 

6j  ■  (ra  *  ra*i)/a  ri  (Vi  ‘  'a-i3  (ra»i  •  V  3  »  1 

The  value  of  b,  follows  the  preceding  formula  for  r^  =  2r.  -  r, ,  where  r, 

-L  0  D  X  D 

is  the  mean  scaled  radius  of  the  inner  duct  wall.  All  the  preceding  co¬ 
efficients  are  computed  in  the  subroutine  MESH. 

2.  The  Vectors:  SUB,  DIAG,  SUPER,  and  RHS 

In  the  mesh  for  the  boundary  value  problem  described  in  Figure  2,  we 
define  the  maximum  number  of  r  ■  constant  lines  by  JMX  (j  )  and  the 

MIX 

maximum  number  of  x  «  constant  lines  by  IMX  (i  ).  For  I  <  10,  the  flow 
field  is  divided  into  a  duct  flow  and.  the  outside  flow  over  the  nacelle 
which  are  solved  separately  in  the  relaxation  process.  The  number  of  points 
along  constant  x.  inside  the  duct  is  JM  (j  in  the  text). 

The  vectors,  SUB,  DIAG,  and  SUPER  are  N  dimensional,  where  N  is  the 
number  of  points  in  a  column.  The  vector  components  are  identified  in  such 
a  way  that 

SUPER(N)  -  SUB(l)  «  0; 
th 

Thus,  the  J  equation  for  each  I  is  seen  to  be  described  by  a  FORTRAN 
statement  of  the  form 

SUB(J)*P(I,J-l)  +  DIAG(J)*P(I,J)  +  SUPER( J)*P(I,  J+l)  -  RHS(j)  (82) 

where  P(I, J)  *  ^  j  and  RHS  denotes  the  right  band  side  terms  of  the 

equations.  The  subroutine  TRISOL  which  solves  the  tridiagonal  system  then 
is  called  by  a  statement  of  the  form 


CALL  TRISOL ( SUB, DIAG, SUPER, RHS,PHI,AUX,N) 


(83) 


The  quantity  PHI  is  the  vector  solution  of  P(I,J)  for  fixed  I  in  the  set  of 
equations  described  by  Equation  (82),  AUX  is  an  M  vector  of  storage  required 
by  the  program.  By  comparing  Equation  (82)  vith  the  equations  in  the 
analysis  one  can  write  down  the  formulae  for  the  vector  components  of  SUB, 
DIAG,  SUPER,  and  RHS. 

Die  quantity  10,  (iQ  in  the  text)  denotes  the  line  x  »  constant  of  the 
mesh  which  intersects  the  nacelle  Una  of  Figure  2  nearest  the  trailing 
edge.  For  all  values  of  i  ^  iQ,  the  column  set  of  equations  is  divided 
into  two  parts,  the  duct  flow  and  outside  stream;  and  the  two  sections  of 
the  column  are  solved  separately  in  the  relaxation  process.  For  i  >  i  , 
the  entire  column  encompassing  both  flows  is  solved  simultaneously. 

3.  Formulation  of  the  Column  Equations 

Because  of  the  initial  condition 

-  f(r)  (84) 

at  x  »  0,  the  column  i  *  1  is  a  special  case.  From  comparison  of  Equation 
(82)  with  Equation  (60),  we  have 

SUB(J)  -  Sj  i  -  2  to  JB 

SUPER(J)  -  b,  J  -  1  to  J  -  1 

J  “  (85) 

DIAG(J)  -  -(v^  +  Sj  +  bj)  i  -.2  to  jm-l 

RHS(J)  ■  "vij  ('i<?aj  -  1  -2  to  V1 

SUPER(JM)  -  0 
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where  f^  *  f(r^)  and 

vu  s  ^**W-*J)  -tr  + 1)  I cu(  q>2i-  9^)  ♦ 

dnfjl 


(86) 


By  comparison  of  Equation  (82)  with  Equations  (45)  and  (49)  we 
obtain 

DIAG(l)  -  -(v11c1  +  bL) 

RHS(l)  -  f  21  -  d1^j_)  + 

DIAG(JM)  -  -(▼ljBC1  +  aja) 

RHS(JM)  -  -V!  (^tf  -  d.fJa)  -  V2R’(x1) 
where  ■*  ^  ^ ^  ") 


(87) 


(88) 


The  J54-1  values  of  are  solved  by  TRISOL  and  iterated.  For  i  =  1 

and  j  j  <.  wt  obtain  similar  relations.  For  J  ■  1  to  JMXL  -  JM  » 

JMX-l-JM,  we  define  the  vectors  for  the  column  above  the  nacelle  as 


SUB( J)  ■  j  SUB(l)  »  0 

SUPER  (J )  -  j 

DIAO(J).  -  (uijBfJ  cx  *  aj^-j  >  b^) 


(89) 


and 


SUPER(JMXl-JM)  «  0 

DUG(l)  .  -  <uljB|i  ex  .  b^)  (90) 

*■<«  ’  “lY  Start  -  Vj.rt>  *  ^  <‘l> 
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where  IL{  j  *  K  —  ( ¥+  t  )  C  c„  l  j  *“  V,  j  ") 


°(  3  "  'jwl/a^aH  ^rJa«-2  "  rJ») 

Since  only  subsonic  inlet  conditions  will  be  prescribed,  the  case 
and  <  0  is  not  considered. 

4.  Formulation  of  the  Column  Equations  2  &  i  ^ 

For  1  <  i  ^  iQ  and  inside  the  duct,  we  have  from  the  comparison  of 
Equation  (82)  with  Equations  (39),  (45),  and  (49),  when  >  0, 


SUB(J)  *  Sj 
SUB(l)  =■  0 
SUFER(J)  -  b 

J 

SUPER(JM)  -  0 


i  -  2  to  A 


J  -  1  to  JB  -  1 


DIAS(JT)  «  -(v^  et  ♦  Sj  ♦  bj)  j  -  2  to  4n  -1 

«®<j)  -  -  VA  .  a  *  -  u>  1  * 5  “  4  • 1 

DIAG(l)  -  -(vu  et  +  bx) 

RBS(l)  -  -^ll(c1^l  +  Xjl  ♦  .u>  ♦*!*'  (*t) 

DIAG(JM)  -  -  (v1J  et  +  a^  ) 

Q  S 

RHS(JM)  -  -  (c^*  +  y  +  difi  -lij  '**2  flL  ^Xi* 

tt  tt  ® 


By  comparison  of  Equation  (82)  with  Equations  (41),  (46),  and  (50),  we 
have  for  the  supersonic  flow,  v^  <  0  and  v^  j  <  0,  the  following  terms 
in  place  of  the  ones  above: 

D1AG(J)  »  -  a^  -bj 

888  ^  *  VlJ  (*i  l^i-lj  *di-l^i-2^;2<J  <‘,m*1 
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APPENDIX  II:  INSTRUCTIONS  FOR  USE  OF  THE  PROGRAM  TRANSDUCT 
1.  The  Input  Data  Cards 

The  parameters  required  by  the  program  are  read  in  by  a  namelist  called 
PARAM.  These  parameters  are  designated  by  the  notation  listed  and  defined 
below: 


Namelist  PARAM 


LSERES 

LRUN 

XMA 

XMJ 


NMAX 

JMX 

IMX 

10 


JM 

MODIN 


Series  number  of  the  flow  configuration. 

Run  number  of  the  series. 

Mach  number  of  undisturbed  exterior  stream. 

Mach  number  of  expanded  fan  jet  flow  far  downstream.  When 
starting  from  an  initial  parallel  flow  (MODIN  *  1),  XMJ 
must  be  less  than  one  for  the  solution  to  converge.  If  a 
supersonic  value  of  XMJ  is  required,  then  a  sequence  of 
solutions  with  XMJ  incremented  by  .05  or  .1  must  be  run 
until  the  value  of  XMJ  is  reached. 

Maximum  number  of  iterations  to  be  computed.  This  is 
chosen  on  the  basis  of  maximum  computing  time  allowed  for 
each  run  and  is  discussed  at  the  end  of  this  section. 

Total  number  of  radial  grid  points  defining  the  mesh  (j 
index) . JMX  -  4l  for  the  examples  described  in  the  report. 

Total  number  of  axial  grid  points  defining  the  mesh  (i 
index).  IMX  ■  60  for  the  examples  described  in  the  report. 

Index  for  the  value  of  x  nearest  the  trailing  edge  of  the 
fan  cowl.  Grid  must  be  constructed  so  that  trailing  edge 
lies  at  midpoint  between  grid  points  X(I0)  and  X(I0  +  1). 

10  is  53  for  the  computed  examples  in  the  report. 

Index  of  largest  radial  mesh  variable  in  the  duct  flow.  JM  - 
20  in  the  example  flows. 

If  MODIN  -  1,  potential  field  is  initialized  for  uniform 
flow.  If  MODIN  ■  2,  potential  field  from  cards  punched  on  a 
previous  run  is  read  for  the  initial  potential  field,  P(I,J). 


MQDCUT 


If  MODOUT  »  2,  the- potential  field  from  the  solution  is 
punched  on  cards.  If  MODOUT  -  1,  no  data  is  punched. 
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Namelist 

Linear 

W1 

CPNAC 


NA 

pxm 


DRFAC 


ERR 
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PARAM 

If  TRUE  (T),  the  linearized  subsonic  flow  is  computed. 

If  FALSE  (F),  the  complete  non-linear  transonic  small 
perturbation  equations  are  solved. 

Over-relaxation  parameter  for  subsonic  colunns .  A  value 
of  1.4  is  suggested.  However,  it  may  vary  between  1  and 
1.9,  but  larger  values  than  1.4  may  lead  to  divergence  of 
the  solution. 

Coefficient  of  pressure  on  the  fan  cowl  at  left  boundary 
of  Fig.  2  for  computing 9 which  in  the  present  formulation 
must  be  restricted  to  subsonic  values.  This  number  is 
usually  taken  from  experimental  measurements.  The  remain¬ 
ing  values  along  x  =  0  in  the  outisde  stream  are  computed 
to  vary  like  l/(r]_ +1)3/2.  This  is  the  same  as  the  incompres 
sible  solution  from  a  source  and  sink  at  x  =  +  1.  (see  p. 
486,  Sec.  16-26  of  Milne-Thomson  (10)).  This-formulation 
was  chosen  for  simplicity  and  convenience  and  should  give 
a  fair  distribution  for  most  calculations. 

Number  of  iterations  between  updating  of Rvalue  at  duct 
input.  When  the  number  is  large,  the  solutions  may  later 
start  to  diverge  if  ip  value  is  not  correct.  Suggested 
values  are  20  to  100. 

Starting  inlet  value  of*>  at  the  fan  plane,  in  the  present 
form  of  the  program  the  distribution  of  <p  is  uniform  across 
the  fan  plane.  This  value  is  updated  every  NA  iterations 
when  WPX  >0.  A  first  estimate  is  required  and  a  reason¬ 
able  value  is  obtained  from  Eq.  (27)  using 

m  =  Af/Ae 

with  the  aid  of  Eqs.  (11),  (l6),  and  (17)  to  determine  K^. 
The  quantity  Af  is  the  duct  area  of  the  fan  plane  and  Ae 
is  the  area  of  the  duct  at  the  transverse  plane  through  the 
fan  cowl  trailing  edge. 

Stretching  factor  for  adjacent  mesh  increments  in  radial 
direction.  Used  for  expanding  the  mesh  in  the  far  field 
where  the  flow  does  not  change  rapidly.  DRFAC  must  satisfy 
the  inequality:  1  ir  DRFAC  ‘*2.  Values  of  1.15  end  1.2 
were  used  in  the  examples  in  the  report.  Its  value  is 
determined  by  how  many  points  are  in  the  grid  and  how  far 
mesh  is  to  be  extended. 

Maximum  error  allowed  for  convergence.  If  the  maximum 
difference  between  two  iterated  values  of  9  becomes  less 
than  ERR,  the  computation  is  terminated;  otherwise,  the 
program  is  terminated  after  NMAX  iterations.  Suggested 
value  for  ERR  is  0.0001. 
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Namelist  PARAM 

WPX  Factor  for  updating  %  at  duct  inlet.  0  £  WPX  —  1. 

If  WPX  -  0,  then  old  value  is  retained.  If  WPX  =  1,  then 

nev  value  computed  from  mass  flow  far  downstream  of  duct 
is  used.  WPX  «  1/2  updates  with  the  arithmetic  mean  of 
former  and  newly  computed  values.  Until  one  observes  how 
the  updated  values  change, every  NA  iteration,  the  value 
of  1/2  is  a  safe  value  for  starting. 

NP  Number  of  iterations  between  intermediate  printout. 


A  typical  set  of  PARAM  data  cards  is  presented  below.  Each  card  is 
punched  starting  in  column  2. 


$FARAM 

M0D0UT-1, 

LSERESy?, 

MODOUT-2, 

LHUN-1, 

LINEAR-.  T., 

XMA-0.8, 

LINEAR—. F. , 

XMJ-0.85, 

wi-1.3. 

NMAX-UOO, 

CPNAC— 0.06 

JMX-41, 

NA-_0, 

DOt-60, 

PXIN—1.21, 

10  33, 

DRFAC-1.1S, 

JIMS, 

ERR-0.0001, 

MQDIN-e, 

NP-20, 

MOD IN-1, 

$END 

For  the  parameters  MOD IN,  MCOOUT,  and  LINEAR,  the  last  of  each  pair 
of  cards  read  is  the  information  used  by  the  program.  Also,  more  than 
one  parameter  may  be  punched  on  a  single  card.  A  space  must  be  left 
between  the  comma  and  the  next  parameter  name. 

After  the  PARAM  list  is  read,  then  a  single  card  with  12  integers 
punched  according  for  FORMAT  (1214)  is  read.  These  are  index  values  of 
tbs  x  grid  at  which  intermediate  values  of  <f  and  (fx  printed  out 

every  NP  iteration. 
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The  coordinates  for  the  duct  and  nacelle  boundaries  are  read  by 
the  subroutine  BODY  along  with  the  axial  coordinates  for  the  mesh  points. 

If  MQDIN  in  the  parameter  list  is  set  equal  to  2,  then  the  cards  for 
P(I,J),  the  values  of  at  the  mesh  points  from  a  previous  run  must  be 

included  at  the  end  of  those  data  cards  read  in  by  the  subroutine  BODY. 

In  summary,  the  data  deck  for  the  program  consists  of  the  following 
(see  Figures  12  and  13): 

1)  Title  card:  TRANSDUCT. 

2)  End  of  record  card  (6,7*8  punches  in  column  1). 

3)  PARAM  cards. 

4)  Card  for  the  twelve  intermediate  values  of  the  x  index  for  which  inter¬ 
mediate  values  of  the  potential  t  and  the  derivative  %  are  printed 
every  NP  times.  Card  FORMAT  is  1214. 

3)  Coordinates  of  the  x  grid  points.  Punched  on  cards  according  to  the 
FORMAT  10F8.6. 

6)  Deck  of  boundary  data  cards.  The  make  up  of  this  deck  will  depend  upon 
how  contours  are  defined  and  data  are  read  by  the  subroutine  BODY 

(see  Figure  13). 

7)  For  MODIN-e,  the  grid  values  punched  from  a  previous  run  must  be 

included. 

For  the  computed  examples  described  in  this  report,  the  running  time 
was  0.61  seconds  on  the  CDC  6600  for  each  iteration  with  41  X  60  »  2460 
grid  points.  Since  the  number  of  computer  operations  is  almost  directly 
proportional  to  the  number  of  grid  points,  then  the  running  time  can 
quickly  be  estimated  for  grids  with  more  or  less  points  than  2460. 
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/XWX  &?/&  &s47X  /vnc/& 2? 
X9GTJ4'<&i/S  XZfiV 


ssc&x  & r&Mfs&y 
erx/?<£>S  __ 

C'SXX  X/G. /<&  -sS- 


FORMAT  C  ^OPS. 6 


IUNCK)  ,  &*!)  12.  .  -po*M*T  02  I  4) 


$£N0 


PARAM 


C£fiD  or  RECORD  CAR# 


TRANSDUCT 
(.TITLE  CARO) 


FIGURE  12  ILLUSTRATION  OF  THE  DATA  CARPS  REQUIRED  BY 
THE  PROGRAM  FOR  MODIN  -  2.  FOR  MODIN  «  1, 
P(I,J)  DECK  IS  OMITTED. 
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FORMAT  C  kO  F  S.O 
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FORMAT  (10  F  8.«)  i 
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FIGURE  13  SUGGESTED  FORM  OF  INPUT  DATA  FOR  BODY  SUBROUTINE 
WHEN  COORDINATES  OF  CENTERBODY,  DUCT  OUTER  WALL, 

AND  FAN  COWL  ARE  USED  TO  DESCRIBE  THE  CONFIGURATION- 
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Subroutine  Body 

This  is  the  subroutine  in  which  the  wall  boundary  is  read  and  the  wall 
slopes  computed  at  the  axial  positions  of  the  mesh.  The  X  or  axial 
coordinates  of  the  mesh  points  are  also  read  in  this  subroutine  with  the 
statement : 

READ  (5..1010)  (X(I),  I  =  1,ZMX) 

1010  FORMAT  (1CF8.6) 

Since  the  contours  for  the  centerbody,  outer  walls  of  the  duct,  and  the 
fan  cowl  are  not  defined  in  any  standard  convenient  coordinate  system, 
this  subroutine  at  least  in  part  may  need  to  be  reprogrammed  for  each  new 
geometry.  The  contours  are  fitted  by  a  spline  program,  then  the  values  of 
the  radial  coordinates  and  the  slopes  at  the  mesh  points  X(l)  are  computed. 
When  coordinates  are  used  to  describe  the  contours,  the  set  up  of  the  deck, 
described  in  Figure  13  is  suggested  with  appropriate  read  statements 
included  in  the  BODY  subroutine. 


The  subroutine  SPLINE  used  by  BODY  to  fit  the  coordinates  is  called  by 
the  statement 

CALL  SPLINE  (M0DE,N,X,Y,D,E,W,JJ,X3,YB,YP,YPP) 

Mode  should  be  set  to  zero  at  first  entry  to  set  up  the  coefficients 
for  fitting  the  N  points  (X,Y)  and  greater  than  zero  thereafter  for 
finding  points  on  the  same  curve.  XB  is  the  value  of  the  abscissa  x  at 
which  the  value  of  the  ordinate  Y3  is  desired.  Since  the  slope  boundary 
conditions  are  required  for  the  main  program,  XB  are  the  values  of  the  x 
grid  points  X(l)  for  which  the  contours  are  defined.  YP  and  YPF  are  the 
first  and  second  derivatives  of  the  fitted  curve  at  the  point  X  =  XB. 

E  and  W  are  storage  required  by  t  he  program  and  both  must  be  at  least  3N 


long.  The  indices  JJ  indicate  whether  end  points  of  the  fitted  curve 
satisfy  the  slope  or  curvature  conditions. 

JJ( 1)  =  0  2nd  derivative  at  left  end  point  given  by  3(1). 

JJ(l)  =  1  1st  derivative  at  left  end  point  given  by  3(1). 

JJ(2)  =  0  2nd  derivative  at  right  end  point  given  by  3(2). 

JJ(2)  =  1  1st  derivative  at  right  end  point  given  by  3(2). 

A  thickness  ratio  DELTA  and  mean  dimensionless  radii  PINNER  for  the 
centerbody  and  RF  for  fan  cowl  must  be  computed  in  the  program  from  the 
contour  data.  RINNER  is  defined  as  the  average  of  the  smallest  and  largest 
radii  along  the  centerbody  and  is  the  radius  at  which  the  linearized  bound¬ 
ary  conditions  are  satisfied.  The  radial  grid  is  constructed  so  that 
PINNER  =  R(l)-(R(2)-R(l) )/2.  The  thickness  ratio  DELTA  is  the  difference 
between  maximum  and  minimum  radii  of  the  centerbody  divided  by  its  length. 
For  the  computed  examples,  this  length  was  conveniently  chosen  as  the.  dis¬ 
tance  from  the  duct  fan  plane  in  Figure  2b  to  the  point  of  maximum  body 
radius . 

Similarly,  the  mean  radius  RF  at  which  the  linearized  boundary  conditions 
on  the  fan  cowl,  on  the  outer  duct  wall,  and  on  the  free  jet  boundary  are 
satisfied  is  the  average  between  minimum  duct  radius  and  maximum  cowl 
radius.  The  radial  grid  is  set  up  so  that  the  midpoint  between  R(JM)  and 
P. ( JM+l )  is  the  radius  RF.  With  RF  and  RINNER  computed  in  the  subroutine 
BODY,  the  main  program  computes  the  radial  grid  points.  Evenly  spaced 
points  are  computed  for  the  duct  flow  and  the  increments  between  radial 
positions  of  the  exterior  flow  are  stretched  by  the  factor  DRFAC  defined 
in  the  namelist  PARAM. 
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For  the  boundary  conditions  required  in  the  main  program,  the  BODY 
subroutine  must  compute  the  derivatives  divided  by  the  value  of  DELTA  for 
the  centerbody,  outer  duct  wall,  and  the  fan  cowl.  The  notation  used  by 
the  main  program  is; 

HX(X),  I  »  1  to  IMX.  Slope  of  centerbody  divided  by  DELIA. 

RLX(I),  I  =  1  to  10.  Slope  of  outer  wall  divided  by  DELTA. 

RUX(I),  I  a  1  to  10.  Slope  of  fan  cowl  divided  by  DELTA. 


In  summary,  the  subroutine  BODY  must  be  programmed  to  compute  RHINE R, 
the  mean  radius  of  the  centerbody;  RF,  the  mean  location  of  outer  duct  wall, 
fan  cowl,  and  free  jet  streamline  boundary;  and  the  slope  boundary  data 
RX,  RLX,  and  HUX  defined  above  at  the  x  grid  points.  The  subroutine  must 
also  compute  the  dimensionless  radii  of  the  centerbody,  RSI;  outer  duct 
wall,  RS2;  and  fan  cow,  RS3,  which  are  required  for  the  subroutine  STRMLN. 


3.  Subroutine  STRMLN 


This  subroutine  computes  the  radial  position  of  the  streamlines  at  the 
x  mesh  positions.  Along  each  radial  mesh  line  R(J)  the  streamline  is 


defined  by 


R1(Xi)  -  (R^) 


jkp 


where  ^  is  the  complete  velocity  potential.  In  terms  of  the  perturbation 
potential,  ,  this  becomes,  since  r  ■  *C  r^  and  'Cfr  ■  ^  . 
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for  points  inside  the  duct  and 


fift 
J  ‘ 

Xf 

R/O-'f’/*)  = 


for  paints  outside  the  duct  flow.  The  derivative  at  the  points  ft(  J ) 
is  found  by  fitting  s  polynomial  through  three  adjacent  points  along  a  line 
x  constant  and  differentiating  it  to  find  ^  .  The  x  integration  is  per¬ 
formed  by  the  trapezoidal  rule. 


4.  Description  of  Printout 


After  the  printing  of  the  parameter  list,  the  coordinates  and  slope  of 
the  lower  and  upper  fan  contours  and  of  the  nacelle  are  printed  for  values 
of  x  corresponding  to  the  mesh  points.  The  actual  slopes  are  printed  in 
the  column  labeled  DER  (far  derivative).  The  slopes  divided  by  the  thickness 
parameter  are  used  as  boundary  conditions  in  the  main  program  and  are 
designated  by  RX,  RLX,  and  RUX  for  the  lower  fan  contour,  upper  fan  contour, 
and  nacelle,  respectively.  The  following  is  a  glossary  of  the  notation  used 
in  the  printout. 

P(I,J)  Value  of  perturbation  velocity  potential  at  X  ■  X(I),  R  *  R(J) 
of  the  mesh 


PHSX  Value  of  ^  at  R(JM). 

FHX(I,J)  Value  of  at  X(I),  R(J). 

R  Dimensionless  radial  variable.  (Multiplied  by 

for  printout  of  PHIX). 

PKIX  Boundary  values  of  %  at  x  «  0. 
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At  every  NP  iteration  through  the  mesh  a  partial  printout  is  made  for 


the  purpose  of  monitoring  the  convergence.  At  this  time,  the  number  of 
iterations  required  for  each  column  is  printed  by  two  rows  of  integers  in 
groups  of  ten.  The  first  row  of  integers  corresponds  to  the  column  itera¬ 
tion  in  the  duct  for  I  ^  10  and  the  entire  column  for  I  >•  10.  The 
second  row  corresponds  to  column  iterating  for  the  region  above  the  nacelle. 
The  maximum  error  is  printed  out  next, along  with  the  I,J  point  of  the  mesh 
at  which  it  occurs. 

At  the  twelve  axial  positions  designed  at  ILINE(K)  two  rows  of  at  J-l 
and  JM  are  then  printed  followed  by  two  rows  of  ^  at  J  ■  JM  and  1. 
The  eighth  row  gives  the  single  quantity  &  at  X  *  X(  10),  the  nacelle 

trailing  edge.  Svery  NA  times  the  value  of  at  the  duct  input  is 

calculated,  printed,  and  updated. 

In  the  final  printout  from  the  subroutine  PRTOUT,  the  column  variables 

are 

CPB  Coefficient  of  pressure  on  the  body  based  on  the  exterior 
stream  conditions. 

CPS  Coefficient  of  pressure  on  the  free  streamline  based  on  the 

exterior  stream. 

PHI  Value  of  along  radial  position  J  ■  JM. 

R-RF  Deviation  of  free  streamline  from  the  trailing  edge  radial 
position. 


REV  SYM 


MO-D6-41078 

*ao«  £8 

i 


By  means  of  the  intermediate  printout  every  NP  times,  the  rate  of 
convergence  of  the  solution  can  be  monitored.  Completely  subsonic  flows  with  as 
initially  prescribed  uniform  flow  (  0,  M0DIN=l)  usually  converge  success¬ 

fully.  When  the  jet  Mach  number  is  chosen  greater  than  unity,  starting 
with  y  =  0  leads  to  an  immediate  arithmetic  error  stop  when  the  program 
encounters  a  negative  index  for  the  P(I,J)  variables.  To  find  a  solution 
for  a  supersonic  free  jet,  the  Mach  number  Mj  (XMJ)  must  be  increased  by 
small  increments  such  as  AM  =  0.05  or  0.1,  starting  with  a  subsonic  value 
and  obtaining  partial  convergence  at  each  Mach  nuniber.  If  the  change  in 
Mach  number  is  too  great,  the  solution  may  take  longer  to  converge  or  it 
may  diverge.  Considerable  computing  time  in  obtaining  a  new  solution  can 
be  saved  by  using  the  P(I,J)  grid  data  from  a  solution  for  which  the 
boundaries  are  similar  to  the  new  configuration  and  the  Mach  numbers  are 
close  to  those  for  the  desired  solution.  Solutions  with  large  imbedded 
supersonic  regions  usually  make  poor  starting  solutions  since  they  often 
cause  the  iterations  to  diverge. 

There  are  three  programmed  stops.  STOPS  1  and  2  occur  with  successful 
completion  of  the  computations.  STOP  1  occurs  for  M3D0UT*1  for  which  no 
data  is  punched.  STOP  2  occurs  for  M0D0UT-2  for  which  the  potential  field 
is  punched  on  cards  as  a  starting  solution  for  a  later  run.  STOP  3  occurs 
when  the  maximum  value  of  the  difference  between  consecutive  values  of 
exceeds  5,  indicating  that  the  solution  is  diverging. 

A  quantity  important  to  the  convergence  is  the  value  of  K,  « 
the  fan  plane  (PXHI  in  PARAM  list).  When  only  a  fair  estimate  is  provided 
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and  the  value  is  not  corrected  in  the  iteration  process,  it  was  found  that 
the  solution  converged  to  a  minimum  value  for  the  error  (ERR)  which  was 
greater  than  0.0001  and  then  started  to  diverge.  Improving  the  value  of 
PXHf  caused  the  solution  to  converge  to  a  smaller  minimum  value  of  ERR 
before  again  starting  to  diverge.  In  the  program,  PXIN  is  updated  every 
NA  times  by  estimating  the  diameter  of  the  jet  far  downstream,  computing 
the  mass  flow  rate  through  this  cross  section,  and  determining  (PXIN) 

which  yields  this  value  of  the  mass  flow  rate  at  the  fan  plane.  The  value 
of  ^  determined  this  way  is  close  to  that  value  of  ^  which  leads  to 
the  most  converged  solution  but  not  necessarily  equal  to  it,  because  of  the 
approximate  way  of  computing  the  downstream  duct  mass  flow. 
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